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We develop a mathematical model for a three-phase free boundary problem in one dimen- 
Oh' si° n that involves the interactions between gas, water and ice. The dynamics are driven by 

melting of the ice layer, while the pressurized gas also dissolves within the meltwater. The 
model incorporates a Stefan condition at the water-ice interface along with Henry's law 
for dissolution of gas at the gas-water interface. We employ a quasi-steady approximation 
for the phase temperatures and then derive a series solution for the interface positions. 
A non-standard feature of the model is an integral free boundary condition that arises 
from mass conservation owing to changes in gas density at the gas-water interface, which 
makes the problem non-self-adjoint. We derive a two-scale asymptotic series solution for 
the dissolved gas concentration, which because of the non-self-adjointness gives rise to 
a Fourier series expansion in eigenfunctions that do not satisfy the usual orthogonality 
conditions. Numerical simulations of the original governing equations are used to validate 
the series approximations. 



Key Words: Free boundaries; Stefan problem; gas dissolution; asymptotic analysis; multiscale; 
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, 1 Introduction 

This paper is concerned with a three-phase free boundary problem involving interactions 
between ice, liquid water, and air. The water-ice interface is driven by a melting process, 
while the gas-water interface is governed by dissolution of gas within the water phase. 
The primary phenomenon we are interested in capturing is the compression or expansion 
of gas that occurs in response to the motion of phase interfaces. 

Free boundaries arise naturally in the study of phase change problems and have been 
the subject of extensive study in the applied mathematics literature [SJ [5J O E] • Mathe- 
matical models of free boundaries generally involve solving partial differential equations 
on some region(s), along with given boundary conditions on a portion of the boundary; 
however, part of the domain boundary remains unknown, and thus some additional re- 
lationship must be provided to determine the free boundary. A classical example is the 

Latest Revision: 1.5 (4 January 2013). Printed: 4 January 2013. 



2 



M. Ceseri and J. M. Stockie 



Stefan problem for a solid-liquid interface [5] that describes a melting or solidification 
process. Here, the primary variable (temperature) is governed by the diffusion equation, 
while the speed of the solid-liquid interface is related to the difference in heat flux on 
either side, which is a statement of conservation of energy. Friedman [6] established well- 
posedness and regularity results for this melting problem, while Crank [5] and Carslaw 
and Jaeger [2j derived analytical solutions using Neumann's method for a variety of phys- 
ical applications. A characteristic feature of all of these solutions is that the speed of the 
free boundary between the phases is proportional to i 1 / 2 , where t is elapsed time. 

Another class of free boundary problems occurs in the study of dissolution and cavi- 
tation of gas bubbles immersed in fluid . IS] . Friedman [7] studied the interface evolution 
for a spherically-symmetric gas bubble immersed in a water-filled container of infinite 
extent, and he proved existence, uniqueness and regularity of the solution. Keller [14] 
studied a similar problem and determined the conditions under which multiple gas bub- 
bles are stable. In particular, he found that gas bubbles should either collapse or else grow 
indefinitely in an infinite medium. In a closed container, however, bubbles can reach a 
stable equilibrium state and he proved that the only stable solution is the one with a 
single bubble. 

This paper was originally motivated by a recent modelling study of sap exudation in 
sugar maple trees during the spring thaw [3|. This is the process whereby maple (and 
other related tree species) generate positive stem pressure that can cause sap to seep 
out of any hole bored in the tree trunk. In late winter there are no leaves to drive 
transpiration, and the maple tree's internal pressure generation mechanism is believed to 
derive from thawing of frozen sap within libriform fiber cells located in the sapwood or 
xylem |17) . These fibers are typically filled with gas during the growing season, but during 
the onset of winter, ice is believed to form on the inner fiber walls thereby compressing 
the gas trapped within. During the spring thaw, the ice layer melts thereby freeing 
the compressed gas which is then free to re-pressurise the xylem sap. A mathematical 
model for sap exudation has recently been developed in the paper [3], which contains 
more details about the physical processes involved. The model predicts build-up of stem 
pressures sufficient to dissolve gas bubbles in the xylem sap, which may also be related 
to the phenomenon of winter embolism recovery that occurs in a much wider range of 
tree species [H] [22] . 

In this paper, we consider a mathematical model for a simpler situation in which a 
closed container is divided into three compartments containing gas, water and ice, in 
that order. While this scenario is not identical to that seen in maple xylem cells, it is 
nonetheless close enough that it permits us to study in detail the dynamics of the free 
boundaries. To our knowledge there has been no other similar study of three-phase flow 
involving gas dissolution and ice melting. There are several other problems arising in 
porous media flow that have some of the same features as our model. For example, the 
modeling of marine gas hydrates [2TJ [24] involves the interplay between gaseous and 
solid hydrates, water, and possibly other components flowing within porous sediments. 
Although these models involve a Stefan condition for a melting front, the gas dynamics 
are driven by hydrate dissociation instead of gas dissolution. Another related problem 
arises in the freezing and thawing of soils contaminated by non-aqueous phase liquids (or 
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Figure 1. Diagram depicting the cylindrical geometry and moving phase interfaces. 

NAPLs) [I3J US] - Here, there is a dissolved gas component but the problem is complicated 
further by the presence of additional phases as well as effects such as mixed wettability. 

The purpose of the present work is to analyze a simple three-phase model that incor- 
porates the dynamics of melting and dissolution. The model is introduced in Section [2] 
and reduced to non-dimensional form. A numerical algorithm is described in Section [3l 
and simulations in Section [4] yield insight into the behaviour of the solution. Motivated 
by these results, we then derive an asymptotic solution in Section [5] that captures the 
essential dynamics, and comparisons are drawn with the full numerical solution. Our 
main aim in this work is three-fold: 

• To understand the basic phase interface dynamics and identify the relevant dimension- 
less quantities and time scales; 

• To develop approximate analytical solutions that can be used either to design more 
efficient numerical schemes or to up-scale material coefficients for microscale models 
such as [3]; 

• To draw connections with existing results on bubble dissolution dynamics. 



2 Mathematical Model 

Consider a cylindrical container of constant radius r and length L (both measured in m) 
that is separated into three compartments containing gas, water and ice as pictured in 
Figure [TJ Assume that the cylinder is long and thin so that L r and we can restrict 
ourselves to a one-dimensional setting where the axial coordinate x varies from to L. 
There are two moving interfaces at locations x — s gw (t) and s W i(t) that separate gas 
from water and water from ice respectively. 

For the sake of simplicity, we consider melting that is driven by a heat source applied 
on the left-hand boundary; and although we will not consider the freezing process, our 
model can be easily extended to handle the freezing case. We are thus interested in the 
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following physical phenomena: (1) heat transfer occuring within and between the three 
phases; (2) phase change at the water-ice interface as ice melts to form liquid water; and 
(3) dissolution of pressurised of gas at the gas-water interface, with subsequent diffusion 
of dissolved gas in the water compartment. The moving boundaries are driven by two 
different mechanisms. The water-ice interface is driven by phase change and the speed 
of the interface is proportional to the difference between heat flux from the adjacent 
compartments (which is the classical Stefan condition [5]). On the other hand, the gas- 
water interface moves in response to changes in volume not only from the dissolution of 
gas in water, but also from the volume change owing to the density difference between 
water and ice. 

We next list a number of simplifying assumptions: 

(Al) The lateral surface of the cylindrical domain is thermally insulated so that heat 
flows only in the axial (x) direction. 

(A2) The system is closed so that the total mass of gas (free plus dissolved) is constant. 
The total mass of liquid and frozen water also remains constant, although the mass 
of the individual phases may change in time. 

(A3) The water and ice densities are constant and are not affected by changes in tem- 
perature. 

(A4) Diffusion in the gas compartment is fast enough in relation to other processes that 
the gas density can be taken as a function of time only. Indeed, considering the 
self-diffusion coefficient for air (D = 2 x 10~ 5 m 2 /s) and a typical length scale 
(d = 100 /jm), the time scale for air diffusion is roughly t d 2 /AD — 1.25 x 10~ 4 s 
(see H sect. 3.32]). 

(A5) The amount of gas that dissolves in the water compartment is small enough relative 
to the initial gas volume that gas dissolution does not significantly affect the motion 
of the gas-water interface. When combined with the previous assumption, this 
implies that the motion of the gas- water interface is due only to the melting of ice 
and the subsequent density increase as ice changes phase from solid to liquid. 

(A6) Neither gas nor water dissolve in or otherwise penetrate the ice layer. 

(A7) Some water is always present in the fiber. Instead of having initial conditions where 
all water is in the frozen state initially, a very thin layer of liquid is assumed to 
separate the gas from the ice. 

We note that many features of this problem are similar to the sap exudation model 
derived in [3] for a closed system that consists of two distinct classes of xylem cells: 
libriform fibers, in which thawing of ice allows compressed gas to force the melted sap 
through the porous fiber wall; and the neighbouring vessels that contain gas and liquid 
sap, where the sap in turn contains both dissolved gas and sucrose. That sap exudation 
model differs from the model we develop here in several respects: 

• we treat only a single compartment containing three phases; 

• in [3], the ice in the fiber is sandwiched between the gas and liquid compartments; 

• we do not consider osmotic effects or sap flow through the permeable fiber/vessel 
boundary; and 
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• the simple cylindrical geometry permits us to neglect surface tension effects due to 
curvature of the gas- water interface. 

In the next three sections, we derive the governing equations and boundary conditions 
in each of the gas, liquid and ice compartments. Following that, we summarise in a 
separate section the remaining interfacial matching conditions that connect solutions on 
cither side of the moving boundaries. 



2.1 Gas compartment 

Denote the temperature in the gas compartment by T g (x, t) [°K] which obeys the diffusion 
equation 

dT q d / dT q \ . . 

for < x < s gw (t) and t > 0. The constant parameters appearing in this equation 
are the thermal conductivity k g [W/m°K] and specific heat c g [J/kg°K], while the gas 
density p g [kg/m 3 ] depends on time according to Assumption lA4[ The initial temperature 
distribution is given by 

T g (x,0) =Tgo(x), (2.16) 

for < x < s gw (0). The heat source that drives the melting of the ice compartment is 
located at the left-hand boundary x = where we impose a constant temperature 

T g (0,t) =Ti >T C , (2.1c) 

that is strictly greater than the melting temperature of ice, T c = 273.15 °K. 

Returning to the gas density, we will now use Assumptions IA2I and IA4I to derive 
a closed- form expression for p g (t) using a mass balance argument that considers the 
total mass of gas (which must be constant) and its division between the gas and water 
compartments. First, the mass [kg] of dissolved gas in the water compartment is 

rs wi (t) 

m w (t) — AM g j C(x,t)dx, (2.2) 

Js gw (t) 

where M g is the molar mass of air [kg/mol], A = irr 2 is the cross-sectional area of the 
cylindrical domain [m 2 ] and C(x, t) is the concentration of dissolved gas [mol/m 3 ] (whose 
governing equation will be given in the next section). The air density may then be written 

as 

_ Ap g (0)s gw (0) - m w (t) + m m (0) 

which along with (I2.2j) determines p g (t) once the dissolved gas concentration and gas- 
water/water- ice interface positions are known. 
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2.2 Water compartment 

We next turn to the water compartment where the temperature T w (x,t) also satisfies 
the heat equation 

(dr w dT w \ a ( 8T W \ 

p wCw ^— +V —j=-i [ k w —), (2.4 a) 

for s gw (t) < x < s w i(t) and t > 0, where c w , p w and k w are the water specific heat, density 
and thermal conductivity respectively. The extra heat convection term on the left hand 
side arises from the slow flow of water due to the melting of ice and the density difference 
between water and ice [SJUHIISS]; the convection velocity v [m/s] will be specified later 
in Section l2~4l We also need to specify an initial temperature distribution 



T w (x,0) = T w0 (x). (2.4 6) 
The dissolved gas concentration C(x,t) obeys the diffusion equation 

oc d ( dc\ 

D w —\, (2.5 a) 



dt dx \ dx r 

where D w is the diffusion coefficient of air in water [m 2 /s]. At the gas- water interface, 
we impose Henry's law 

C(s gw (t),t) = ^ Pg (t), (2.5 6) 

which states that the concentration of gas dissolved at the interface is proportional to the 
density of the gas in contact with the liquid. Here, H denotes the dimensionless Henry's 
constant. Finally, we impose the initial condition 

C(x,0) = C (x), (2.5 c) 

and the following no-flux boundary condition at the water-ice interface 

dC 

— ( Swi (t),t)=0, (2.5 d) 

which is a simple statement of the fact that dissolved gas does not penetrate the ice (in 
accordance with Assumption IA6|) . 

The following two remarks relate to the distribution of air between the gaseous and 
dissolved phases. 

Remark 1 (Conservation of air) We first show that (|2.3[) implies conservation of mass 
for total air in the gaseous and dissolved phases. The total mass of air at any time t can 
be written as the sum of the air in the water and gas compartments: 

m — m w (t) + A / p g (t)dx, 
Jo 

= m w (t) + As gw (t)p g (t). 
Then, replacing p g {t) with (|2.3p leads to 

m = m. w (Q) + Ap g (0)s gw (0). 
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This last expression is simply the sum of the total initial mass of dissolved and gaseous 
air, and hence the total mass m of air is conserved. 

Remark 2 (Connection with Keller's analysis of gas bubble dynamics) Our aim here is 
to derive an expression for the rate of change of the gas density, which can then be related 
directly to an equation derived by Keller for the dynamics of dissolving gas bubbles in 
water /14/. To this end, we take the time derivative of the gas density from equation 

dp g -m w (t)s gw (t) - s gw (t) [Ap g (0)s gw (0) - m, w (t) + m w (0)] 
~dt ~ As*(t) ' 



m w (t) Pg(t)Sg W {t) 



ASg W (t) 



,(t) 



(2.6) 



m w (t) — AM g 



m w {t) = AMg s wi {t)C{s wl {t),t) - s gw (t)C{s gw (t),t) + D w {^-(s wl i! ). / ) 



where the "dot" denotes the time derivative. An expression for rh w (t) can be obtained by 
differentiating (12.21) 

rs mi (t) QQ 

s W i{t)C(s wi (t),t) - s gw (t)C(s gw (t),t) + / — (x,t)dx 

Js gw (t) ot 

The integral term containing d t C may be integrated directly by first replacing dtC using 
the concentration equation (12.5 a\f 

'dC dC 

dx 

and then applying the boundary condition (|2.5 



(Sg W (t),t) 



r dC 
m w (t) = AMg s wi {t)C{s wi (t),t) ~ s gw (t)C(s gw {t),t) - D w — (s gw (t),t) 



We may now substitute this last expression for rh w (t) into equation 



to obtain 



1 d(p g s gw ) 



M„ 



dt 



dC 

■jC(s gw , £) s w iC {s w i , t*) -\- D w \S g w ? (^*^) 



This equation coincides with Keller's equation (2.6) /14j /. except for slight differences 
arising from to the fact that we are working in a cylindrical geometry and we also include 
the water-ice interface motion in the evolution of the gas density. 



2.3 Ice compartment 

In the ice compartment, the equation governing the ice temperature Ti(x,t) is 

dT, d (. dT t 



dt dx \ dx J ' 



(2.8 a) 



for s W i(t) < x < L and t > 0, where Cj, pi and fcj are the specific heat, density and 
thermal conductivity of ice. The initial temperature distribution is given 

r i (aj,0) = T i0 (x) ) (2.8 6) 
and on the right boundary we impose a convective condition of the form 

8T 

-k i -±(L,t) = e(T i -T 2 ), (2.8 c) 
ox 
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where 9 is a convective heat transfer coefficient [W / m 2 °K] and T2 is a given ambient 
temperature. 



2.4 Interfacial and matching conditions 

We now state the matching conditions at the two phase interfaces. First, we require that 
the temperature and heat flux are both continuous at the gas-water interface 

T g {s gw {t),t) = T w (s gw (t),t), (2.9) 

dx ^ Sgw ^ ' ~ dx 
Based on geometric and conservation arguments, we can derive an equation for the evo- 



■{s gw {t),t) = k w -^{s gw {t),t). (2.10) 
v 

lution of the gas- water interface 



s gw (t)= (l-y)s W i(t), (2.11) 

which relates the velocities of the two interfaces via the difference in volume owing to 
contraction and expansion of ice (details of the derivation are provided in Appendix |A"]) . 
This condition is also given in Crank's book [51 Eq. 1.32] which expresses the velocity of 
the liquid phase when a density change is taken into account; therefore, we impose 

v = s gw (2.12) 

for the convection speed in equation (|2.4 al) . 

At the water-ice interface, the temperature must be continuous 

T w (s wi {t),t) = Ti{s wi {t),t) = T c , (2.13) 

with the added requirement that the temperature on both sides of the interface must 
equal the melting point. The evolution of the water-ice interface is governed by 

dTi dT w 

^PiSwi — k{— k w — — at x = s w i(t) 7 (2.14) 

ox ox 

where A is the latent heat of melting per unit mass [J/kg]. This Stefan condition is a 
statement of conservation of energy, where the amount of heat generated by the change 
of phase (\piS W i) is balanced by the difference in heat flux from either side of the phase 
interface. Finally, to close the system we require initial conditions for the phase interface 
locations: 

Sgio(O) = s gw0 , (2-15) 
s wt (0) = s wi0 . (2-16) 



2.5 Parameter values 

The values of all parameters defined above are given in Table [T] in SI units and are 
taken from the data for the sap exudation model in .3;. The geometrical parameters 
L = 10~ 3 m, r = 3.5 x 10~ 6 m and A = -rrr 2 = 3.85 x 10 -11 m 2 are all based on the size 
of a typical libriform fiber in the xylem of a maple tree. 
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Parameter 


Symbol 


Units 


Value 


Domain length 


L 


m 


1.0 x 10" 3 


Domain radius 


r 


m 


3.5 x 10-6 


Cross-sectional area 


A = 7T7- 2 


m 2 


3.85 x lO^ 11 


Densities 


Pg 


kg/m 3 


1.29 




Pw 


kg/m 3 


1000 




Pi 


kg/m 3 


916 


Specific heats 


C 9 


J/kg°K 


1005 




c w 


J/kg°K 


4180 




Ci 


J/kg°K 


2050 


Thermal conductivities 


kg 


W/m°K 


0.0243 






W/m°K 


0.58 




k-i 


W/m°K 
m 2 /s 


2.22 


Thermal diffusivities, a = k/(pc) 


OCg 


1.87 x 10~ 5 






m 2 /s 


1.39 x 10~ 7 




OLi 


m 2 /s 


1.18 x 10~ 6 


Diffusivity of dissolved air in water 


D w 


m 2 /s 


2.22 x 10~ 9 


Molar mass of air 


M g 


kg/mol 


0.0290 


Henry's constant 


H 




0.0274 


Convective heat transfer coefficient 





W/m 2 °K 


10.0 


Latent heat of melting 


A 


J/kg 


3.34 x 10 5 


Critical (melting) temperature of ice 


T c 


°K 


273.15 


Left boundary temperature 


Ti 


°K 


T c + 0.005 


Right boundary temperature 


T> 


°K 





2.6 Non-dimensional equations 

We now non-dimensionalise the governing equations by introducing the following dimen- 
sionless quantities denoted by a superscript asterisk (*): 

x Lx , Sg W = LSg W , s w i = Ls wi , 

t = tt*, C = CC*, p g =p g p* g , (2.17) 

T^Tc + i^-TJT*. 

We choose as a length scale L = 1CP 3 m, which corresponds to the typical length of a 
libriform fiber [3J. Density is rescaled by the value p g = 1.29 kg/m 3 for air at 1 atm 
and 0°C, and the concentration by C — p g /M g . The time scale i is chosen equal to the 
characteristic scale typical in Stefan problems for motion of the water-ice interface 

t= , *\ y (2.18) 
k w (Ti - T c ) 

because the melting process is the driving mechanism for this problem. 

Substituting the expressions from (|2.17|) into the model equations and dropping as- 
terisks to simplify notation, we obtain the following dimensionless system where all new 
parameters are listed in Table [21 
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Table 2. Characteristic scales and dimensionless parameters. 



Parameter Expression Units Value 

Pa kg/m 3 1.29 

C mol/m 3 44.6 
M a 

L 2 

t w — s 7.21 



to 



OLw 



T 2 

U — s 8.46 x 1CT 1 

cm 

- L 2 \p l t w 6 5 

s 1.06 x 10 



k w (Ti-T c ) St 



St 

f 2 



Pi 



T 2 -T c 
Ti-T c 

Olgt CXg [3 yj 

L? a w 
a w t 8 
~L? = St 



0.916 



V ^ 23.9 

Kg 

V 3.83 

Bi — 4.50 x 10" 3 

Le ^ 62.5 



(Tl ^ c)Cw 6.26 x 10" 5 

A 



1.97 x 10 6 
1.46 x 10 4 



% = ^ 1.25 x 10 5 



In the gas compartment, < x < s gw (t): 

T ff (x,0) =T s0 (a;) for < a; < * flU) (0), (2.19 6) 

T fl (0,t) = l, (2.19 c) 

where the dimensionless diffusion coefficient (3 g = a g t/L 2 and a g = k g /(p g c g ) is the 
thermal diffusivity of air. 

In the water compartment, s gw (t) < x < s w i(t), we have equations for both tempera- 
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ture and concentration 

BT BT B 2 T 

u± w , u± w u ± w . . 

T w (x, 0) = T w0 (x) for s gw {0) <x< s m (0), (2.20 b) 



dC 6 d 2 C 



(2.21a) 



dt St Le Bx 2 ' 

C{x,0)=C (x) ior s gw (Q) <x <s m (0), (2.216) 

C(s gw (t),t)=Hp g (t), (2.21c) 
8C 

— (s m (t),t) =0. (2.21 d) 

Here, f3 w = a w i/L 2 , St = c w (T\ — T c )/\ is the Stefan number and the Lewis number 

Le = a w /D w is a dimensionless ratio of thermal diffusivity of water to the diffusivity of 
dissolved gas. 

In the ice compartment, s W i(t) < x < 1, 

dTi 3 2 T t , 

af = A7^, (2 ' 22a) 

ri(a?, 0) = T M {x) for s wi (0) < x < 1, (2.22 6) 

9T- ~ 

-^(l,t) = Bi(T,(l,t)-T 2 ), (2.22 c) 

where ft = aJ/L 2 and T 2 = (T 2 - T c )/(Ti - T c ). The Biot number Bi = LO/h is a 

measure of the relative resistance to heat transfer of the outer surface of the ice to that 
in the interior. 

The non-dimensional forms of the intcrfacial conditions at the gas-water interface are 

BT BT 

^(s gw (t),t)=V-^r(s 9W (t),t), (2.23) 

T g (s gw (t),t) = T w (s gw (t),t), (2.24) 
while at the water-ice interface 

T w {s wl (t),t) = Ti{s wi (t),t) = 0. (2.25) 
Equation (|2 .3[) for the gas density reduces to 

p g {t) = J -^B Js - (t) . (2.26) 

s gw\t) 

The gas-water interface equation (12.11[) becomes 

Sgw(,ty (1 $wi> 
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where S = pi/p w , which can be integrated directly to obtain 

s gw (t) = Ai+ A 2 s wi (t) \ 
where A 1 = s gw {0) -(1-6) s wi (0) \ . (2.27) 
and A 2 = 1 — S J 

Finally, the dimensionless form of the Stefan condition (|2.14l) is 

Swi = \^^^{s W i,t) ~ ~dx'^ Sm,t ^) ' (2.28) 

where -0 = ki/k w . Note that the choice of time scale i made in (|2.18[) was made so 
that the coefficient in front of s W i scales to one. Comparing i the typical sizes of the 
corresponding scales for heat diffusion (t g , t w and ti in Table [5]) , it is clear that the front 
motion occurs over a much slower time scale. 

In summary, our model consists of a coupled nonlinear system of equations that is 
composed of: 

• four PDE initial-boundary value problems (|2.19|) - (|2.22[) for the temperatures and dis- 
solved gas concentration; 

• one ODE initial value problem (|2.28j) for the water-ice interface position; 

• two algebraic equations (|2.26[) and (12.27)) for the gas density and gas-water interface 
position. 

Because of the nonlinearity present in the equations, it is not possible to derive an 
explicit analytical solution and so we must resort to numerical simulations or approximate 
analytic methods. In the next section, we describe our numerical discretisation procedure 
and present approximate results that in turn suggest an appropriate choice of analytic 
solution. 



3 Numerical solution algorithm 

We begin by briefly describing our approach for solving the PDEs governing temperature 
and concentration. We use the method of lines, discretising the PDEs in space on a cell- 
centered grid and then solving the resulting system of time-dependent ODEs. In order 
to capture moving boundaries sharply, we employ a moving mesh approach in which N 
equally-spaced grid points are distributed over each of the gas, water and ice domains, 
so that 

X J g(t) = (j - l/2)/l fl (t) With MO 

4u(t) = s gw (t) + (j - l/2)h w {t) with h w (t) 
xl(t) = s wi [t) + {j - l/2)fti(t) with hi(t) 

for j — 1,2, ... ,N, and where h^(t) for t = g,w,i denotes the grid spacing on the 
corresponding compartment. When using such a moving computational grid, we must 
introduce an additional convection term in each parabolic PDE owing to the grid mo- 



s gw jt) 
N 1 

' N ' 

N ' 
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tion [TOl [12] 

df df d 2 f 

~dt~ u fa- K dxt' (3 ' 1} 

where f — T g , T w , Ti, C and k = /3 g , fi w , Pi, <5/(StLe) respectively. The convective term 
has a velocity u that corresponds to the mesh velocity xt on the gas and ice compartments, 
but equals x w — s gw on the water compartment owing to the presence of the convective 
term in (|2.20 a\i . 

The spatial derivatives appearing in equation (|3.1[) are replaced using centered, second- 
order difference approximations to obtain 

~dt u * 2h t ~ Kl hj ' [3 - 2) 

where n(t) f(x^,t) are the discrete approximations of the dependent variables for 
£ = g,w,i and j = 1,2, ... ,N. Centered finite differences are also used to discretise 
the boundary conditions and hence maintain second order accuracy throughout. This 
requires values of the approximate solution that lie within a grid cell lying immediately 
outside of each compartment; to this end we introduce fictitious points } — he and 

x^ +1 — x 1 ^ + hi. A Dirichlet boundary condition such as (|2.19 c\ is approximated using 
an arithmetic average 

9 _ 9 _ i 

2 

which is solved for the fictitious value as T° = 2 — T g . Furthermore, a Neumann boundary 
condition such as (|2.22 c[) is approximated by 

JiiV+l _ rpN 




which yields 



2-h l Bi)T» +2hjBiT 2 
2 + hiBi 



The remaining fictitious point values are obtained in a similar manner using the other 
boundary and matching conditions. Finally, integrals of concentration that appear in 
the boundary condition (|2.21 c|) (via the density (|2.26[1 ) are approximated using the 
trapezoidal rule, so that the resulting spatial discretisation is fully second order in space. 

The semi-discrete temperature and concentration equations comprise a system of AN 
time-dependent ODEs. One additional ODE derives from the water-ice interface equation 
(|2.28[) in which the spatial derivatives are also approximated using centered differences. 
The resulting system of AN + 1 ODEs is implemented in the Matlab® programming 
environment and integrated in time using the stiff solver odel5s. In all cases, the error 
tolerances for odel5s are set to AbsTol=le-10 and RelTol=le-8. 
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4 Numerical simulations 

The method described in the previous section is now employed to simulate the model 
equations and to evaluate its sensitivity to various physical parameters. In all simulations, 
we make the following choices for initial conditions: 

• s gw (0) = 0.1 and s W i(0) = 0.11, so that the water is initially completely frozen except 
for a thin liquid layer (refer to Assumption IA7I) ; 

• C(x, 0) = 0, corresponding to no dissolved gas; 

• T g (x, 0) = T w (x,0) = 1 and Tj(x, 0) = T2, so that the gas and water compartments 
are both equilibrated with the left boundary temperature, while the ice compartment 
is equilibrated with the ambient (sub-freezing) temperature at the right boundary. 

We begin by focusing on the water-ice front motion that drives the phase change 
dynamics and in turn influences the gas dissolution. We consider as a "base case" the 
situation where the left and right boundary temperatures are T\ = T c + 0.005 and 
T2 = T C — 0.005, for which results are provided in Figure [5J Two other cases with larger 
values of T\ and T2 are presented in Figures |3] and HI for comparison purposes. Note that 
all plots are show in dimensionless variables. 

The the base case, the plot in Figured of the dissolved gas concentration (measured 
at the right-hand boundary) exhibits a clear division of the solution behaviour into three 
separate time periods: 

(1) A very short initial transient during which the concentration undergoes a rapid 
increase from zero at t = to some maximum value at t ~ O(10~ 7 ). This tran- 
sition layer arises because we have chosen initial conditions corresponding to zero 
dissolved gas and hence Henry's law forces the initially very thin water layer to 
rapidly "fill up" with gas. The corresponding diffusion of dissolved gas within the 
water compartment is easily seen in Figure [2]:. 

(2) The gas concentration remains roughly unchanged over the interval t £ [10 -7 , 10 -2 ] , 
since the water does not yet melt appreciably. 

(3) The time t m 10~ 2 signals the onset of ice melting, after which the water com- 
partment begins to grow in size. Even though this allows more gas to dissolve in 
the water layer, the increased volume leads to a decrease in the dissolved gas con- 
centration as the pressure in the gas compartment decreases. This effect is evident 
from Figure^! where we observe that the concentration profiles through the water 
layer are roughly constant in x, although there is a very slight increase in C from 
left to right. 

The presence of these three, clearly separated time scales is a characteristic feature of 
the evolution of dissolved gas. Because the concentration profiles are almost constant in 
x over longer times, the gas concentration dynamics are driven primarily by the relative 
motion of the free boundaries. 

A somewhat counter-intuitive result derives from the observation that after initial 
transients are complete, C attains its maximum value at the water-ice interface rather 
than at the gas-water interface where dissolution is actually taking place. This slight 
positive slope in the plot of C versus x becomes more pronounced as the boundary 
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temperature difference T± — Ti is increased, and can be seen most clearly in Figure [4]d 
where the temperature difference is largest. We also remark that the speed of the free 
boundaries increases with T\ — Ti which allows less time for the gas to adjust in the 
water compartment. 

We close our discussion of the base case with a look at the final two plots in Figure [21 
The water-ice interface in Figure [5J1 shows the expected sub-linear behaviour that is 
consistent with the t 1 / 2 dependence predicted by the analytical solution to the Stefan 
problem. This behaviour is confirmed by our asymptotic results in Section[52I According 
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(a) Gas concentration at x — L. 
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Figure 3. Solution plots with boundary temperatures T\ = T c + 0.005 and T 2 = T c — 0.02, 
with T 2 = —4. 



(a) Gas concentration at x — L. 
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(b) Gas concentration profiles (long time). 
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Figure 4. Solution plots with boundary temperatures T\ = T c 
f 2 = -1. 



1 and T 2 = T c - 1, 



to Figured, the temperature field is a continuous function that changes relatively slowly 
over time. Furthermore, the temperature is approximately linear within each compart- 
ment, with a pronounced "kink" at the each interface locations. Both of these results will 
be explained by the analytical solution derived in Section O 

The effect of increasing the temperature difference T\ — T2 can be seen by comparing 
the results in Figures [3] and 0] with Figure [5] There is a significant slowing of the ini- 
tial transient gas dissolution dynamics as T\ — T2 is increased, although the long-time 
concentration dynamics are largely unchanged. However, as mentioned above, there is a 
slight increase in the slope of the concentration profiles in Figure [4]b. 

Comparing Figures [2HIJ we remark that for the base case with temperature increments 
of 0.005, the ice layer melts away after about one hour. In contrast, the melting time 
shortens to 17 seconds when the temperature increment is taken as large as 1.0. The 
only place that T\ enters the model is through the Stefan number St, which explains why 
changes in T\ have the effect of altering the time scale for the free boundary motion. 

We conclude this section by investigating the effect of taking a relatively large initial 
value for the dissolved gas concentration, C(x, 0) = 0.055, rather than taking C(x, 0) = 



Three-phase free boundary problem with melting and dissolution 17 
(a) Concentration profiles (short time). (b) Gas concentration at x — L. 




as we have so far. We will see later on that this initial concentration is large in the sense 
that it is twice the steady-state value of concentration for the base case. Hence, this 
situation may be viewed as corresponding to a "super-saturated" case in which one would 
expect dissolved gas to immediately cavitate and form bubbles. The results in Figure [S£i 
are consistent with this hypothesis, and show that the behaviour of the concentration 
profiles is reversed relative to the base case in Figure [3b, in that concentration decreases 
from the initial value to its quasi-steady state. 



5 Approximate analytical solutions 

Motivated by the numerical results in the previous section, we now derive an approximate 
analytical solution that is based on the following observations: 

• The temperature is approximately linear within each compartment, and equilibrates 
rapidly to any change in conditions over the time scale of the interface motion. This 
suggests using a quasi-steady approximation for each temperature variable. 

• The dissolved gas concentration evolves over two distinct time scales: a rapid initial 
equilibration phase driven by diffusion (on the order of 10 _8 -10 -5 seconds) during 
which gas dissolves at the gas- water interface to fill the liquid compartment; and a 
much longer time scale corresponding to the onset of ice melting (on the order of 
t = 0.01-0.1 s) when the water-ice interface begins to move and the volume of the 
water compartment increases appreciably. 

As a result, we approximate the solution in three stages. First, we make a quasi-steady 
approximation for temperature that permits us to write Ti(x, t) as linear functions of x 
for £ — g,w, i, that vary in time only through changes in the interface locations. Second, 
we derive a simpler ODE for the water-ice interface s W i(t) that makes use of a series 
expansion in the small parameter Bi, which then also yields an approximation for s gw (t) 
via equation (|2.27|) . Finally, we develop a two- layer asymptotic solution for the dissolved 
gas concentration C(x,t) based on the separation of time scales mentioned above. 
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5.1 Quasi-steady approximation for temperatures 

The time scales for diffusion of heat in the gas, water and ice compartments can be 
estimated using 

L 2 

f,£ = — for £ = g,w, i, 
on 

where the thermal diffusivities ai and length scale L are taken from Tables Q] and [2l The 
corresponding time scales are t g ps 5.4 x 10~ 4 s, t w w 7.2 x 10~ 2 s, and s» 8.5 x 10~ 3 s. 
In contrast, the time scales for motion of the gas- water and water-ice interfaces were 
observed in the numerical simulations from the previous section to be at least one order 
of magnitude larger than this; consequently, the phase temperatures will adjust rapidly 
in response to any motion of the interfaces. It is therefore reasonable to assume that the 
temperatures Tg are quasi-steady in the sense that they do not depend explicitly on time 
but instead have an implicit dependence on t through the free boundary locations s gw (t) 
and s m (t). 

The convective term in the water equation (|2.20 a[) is so small (on the order of 10~ 2 ) 
that it is reasonable to neglect. Therefore the temperature equation in all three com- 
partments has the simple form d xx Tg — and consequently the temperature is well- 
approximated by linear functions of x 

T e (x,t) = at(t)x + b e (t) for £ = 5, w, i. (5.1) 

The coefficients ag(t) and bg(t) can be determined by imposing boundary and matching 
conditions (12.19 c|) . (|2.22 cl) and (|233| - (|2~25]) . after which we obtain 

T (x t) = — - Sgw ^ + Ssw - Swi (5 2 a) 

T w (x,t)= °y-_ X , (5.2 6) 

1 + Bi(l - s wi ) 

on the corresponding sub-intervals. 



5.2 Asymptotic expansion for water-ice interface 

We next derive an analytical solution for the water-ice interface s W i(t) by substituting 
the approximations just derived for T and T w into the Stefan condition (|2.28|) along with 
the expression (I2.27|) for s gw to obtain the following ODE 

{B l + B 2 s wi )(l + Bi(l - s wi ))s wi = (1 + Bi(fl 5 + B 6 s wi )). (5.3) 

The constants appearing in this equation are 

B 1 = (r]-1)A 1 , B 2 = l + (r]-l)A 2 , 

B 3 = B lSwi (0) + £ 2 2^M B 4 = B lSwi (0) + 5a=li a a.( ) - ^4,(0), 
B 5 = 1 + il)T 2 B x , B 6 = ^T 2 B 2 - 1, 
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while A\ and A 2 are the same constants defined earlier in equation (|2.27p . This ODE 
can be integrated in time over the interval [0, t] to obtain the following integral equation 
for s wi : 



Bis h 



B 



Bi 



R i B2 ~ Bl 



Bo 



2 

= B 3 



s 2 ■ - — s 3 



i + BiB 4 + Bi 



B 5 t + B e 



i(l)dl 



(5.4) 



Because the Biot number satisfies Bi <C 1 (see Table [2]) it is reasonable to look for a 
series solution of the form 

s wi (t) = s (t) + Bi»i(t) + 0(Bi 2 ). (5.5 a) 

Substituting this expression into (I5.4[) and collecting terms in like powers of Bi, we find 
that to leading order 

1 

B~2 



*o(t) 



B\ + 2B 2 (B 3 +t)-B 1 



(5.5 6) 



while the next order correction is 

1 



«i(t) 



Bo , . q B] — Bo , . q _ _ 

,> , M i^o(t + 9 2 s (t) 2 - B lSo (t) + B 4 

Bi + B 2 so(t) \ 3 2 



+ B 5 t + B 6 / s a (l)dl 



(5.5 c) 



Using the water-ice interface approximation in equations (I5.5[) the gas-water interface 
may be determined from (|2.27p . 

We conclude this section by drawing a connection between the leading order solution 
So{t) in the limit as Bi — > and the classical solution of the Stefan problem where the 
melting front moves with a speed proportional to t 1 / 2 . Although equation (15.5 b\ does 
not have exactly this form, the behaviour is consistent in the limits of large and small 
time. In particular, if we expand (|5.5 b\ in a Taylor series about t = we find that 



so(f) 



^B\ + 2B 2 B 3 - B l 



B 2 ^JB\ + 2B 2 B 3 

Furthermore, the large-time limit of (15.5 61) yields 



t + 0(t 2 ) (as i ->()). 



(5.6) 



so(t) 



-Bi 

B 2 



2t 

B~2 



1/2 



B 2 + 2B 2 B 3 



2B 2 



2t 
B 2 



-1/2 



^- ) +0{t- 3 / 2 ) (aat->oo). (5.7) 



When these two series expansions are plotted against the exact expression for so(t) in 
Figure El we see that both match well for small and large times, and in particular the 
large-time expansion (|5.7|) shows the expected t 1 / 2 behaviour. 



5.3 Two-scale asymptotic solution for gas concentration 

The numerical simulations from Section |4] (more specifically, the plots in Figures [2ji, [3^,, 
Hk) exhibited a clear separation of time scales during the evolution of the dissolved gas 
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Figure 6. Series expansions of So(t) for large and small times, showing the expected y/i 
behaviour as t — > oo (note that the i-axis is on a log scale). 



concentration. Starting from the given initial value, the concentration increases rapidly as 
gas dissolves at the gas- water interface and diffuses throughout the water compartment. 
We repeat our earlier observation that the gas concentration is nearly constant in space, 
but has a slight positive slope that leaves the maximum value at the water-ice interface 
(see Figure 0J}) ; this maximum is achieved over the short diffusion time scale before the 
free boundaries begin to move. From then on, the dissolved gas concentration remains 
essentially linear and decreases over a much longer time scale that is driven by the motion 
of the free boundaries. It is this dual time scale behaviour that we aim to explain in this 
section. 

To this end, it is helpful to derive rough estimates of the time and length scales involved. 
The time required for the dissolved gas to diffuse a distance d — j(s w i(0) — s gw (0)) 
corresponding to half the width of the water compartment is 

tp 

t d = — m 1.13 x 1CT 2 s. (5.8) 

This value should be compared with the time t W i required for the ice to melt completely, 
which can be estimated by setting s W i(t W i) = 1 in equation (|5.5I) and focusing on the 
leading order term to obtain 

twi = r(5 a + 2S i - 2B z) ~ 96 A hours ' 

which is six orders of magnitude larger than the diffusion scale td in (|5.8[) above. Moreover, 
over this same time scale, the water-ice interface is only capable of travelling a distance 
of 

L(s m (t d ) - s wi (0)) ps 2.25 x 10- 8 L. 

Hence, the phase interfaces can certainly be treated as stationary over the diffusive time 
scale td- 
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Based on these observations, we now develop a two-layer asymptotic expansion for the 
dissolved gas concentration. We begin by rescaling the dimensionless spatial variable 
according to 

X- 3^,(0) 

As 

where As = s w ,(0) — s gw (0) > 0. By substituting into equation (|2.21 al) and defining a 
new concentration variable G(y,i) = C(x,t), we obtain 

where the new diffusion parameter is 

LeSt.. . 9 
e = ——(As) 2 < 1. 
o 

It is convenient at this point to rescale the interface positions according to 

/,N S Wi{t) — Sg W (0) Sg W (t) — Sg W (0) 

&wi(t) = - and o- gw [t) = . 

As aw ^ s 

Two series expansions will next be developed for the concentration variable G(y,t): one 
on an "outer region" corresponding to times t = 0(1), and the second on an "inner 
region" corresponding to t = 0(e) <C 1. 



5.3.1 Outer expansion (large time) 

For large times, we suppose that the dissolved gas concentration is a series in the small 
parameter e: 

G(y,t) - G (y,t) + eG^y, t) + 0(e 2 ). (5.11) 

Substituting this expression into (|5.10|) and collecting terms with like powers of e gives 
rise to the leading order equation 

dy 2 

which has solution G(y,t) — a(t)y + b(t), similar to the quasi-steady approximation for 
temperature we obtained in Section 15.11 The leading order boundary conditions corre- 
sponding to (|2.21 c|) and (|2.21 rfj) are 

^—^-(a wi (t),t) = 0, 
dy 

( /■! r^dt) \ 

H( + H[ / C (y)dy- / G (y,t)dy 

O {a gw {t),t) = — — , 

where we have introduced the notation 

C=^, (5.12) 
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which is a positive constant because As > by Assumption IA7I The zero Neumann 
boundary condition requires that a(t) = 0, after which we obtain the leading order 
solution 



G (v,t) 



HQ + h( C (y)dy 
Jo 



(5.13) 



C + °~gw + H(<T wi — Og W ) 

At the next higher order in e, we obtain the following boundary value problem for 
Gi{y,t) 

<9 2 Gi dG 



BGi 
dy 



dy 2 

(a wi (t),t) = 0, 



at ' 



G 1 (a gw (t),t) = -H^(t) 



<T sro (t) 



Gi(y,t)dy, 



where we have defined 



( + (Tg W (t)' 

Using a similar argument to the leading order solution, we obtain 



G 1 (y,t) = ^(y,t) 



y 



(5.14) 



(5.15) 



Note that dtGo < so that G\ is an increasing and concave downward function of 
y that attains its maximum value at the right-hand endpoint y = a w f, therefore, the 
asymptotic solution exhibits the same behaviour observed earlier in the numerical results 
for concentration in Figure HJd. 



5.3.2 Inner expansion (small time) 

For much shorter times with t = 0(e), we rescale the time variable according to 

t 



(5.16) 



and also denote the inner solution for concentration by j{y,r) — C(x,t), where y is 
the same rescaled spatial variable in (|5.9[) . Under this scaling the concentration diffusion 
equation (|5.10|) reduces to 

As mentioned before, over such a short time interval the phase interfaces are essentially 
stationary so that we can look for a solution 7 on the fixed interval y € [a gw (0), a W i(0)] — 
[0, 1]. The initial and boundary conditions ([2.21 6|) - (|2.21 a} may then be written in terms 



Three-phase free boundary problem with melting and dissolution 



23 



of 7 as 

7(V,0) = C (V). 



H / ,1 ,i 



7 (0,T) = ff+- M C (y)dy-^ 1 (y,r)dyj, 
| I (l ! r) = 0. 

We begin by determining the steady state solution for this problem, which is simply 
the constant value 



H( + Hf C (y)dy 
Jo 



C + H 

We then define j(y, r) = 7(2/, r) — 700, which satisfies the same equation f|5.17[) along 
with the following modified initial and boundary conditions 

l(y, 0) = C {y) - 7oo, 

7(0,t) = --r A /{y,r)dy, 
s. Jo 

^(l,r)=0. 

This modified problem can be solved by the method of separation of variables to obtain 

00 

j(y, t) = Y^ a n cos (^„(y - l))e^- T , (5.18) 

n=l 

where fi n are solutions to the nonlinear equation 

H n (, + H tan /x n = 0. (5.19) 

In the method of separation of variables, it is customary to determine the series coef- 
ficients a n by multiplying the initial condition 



Co{y) - 700 = ^2 0,1 cos {^{y - 1)) 



71=1 

by another eigenfunction from the set T = {cos(/i„(y — 1)) | n = 1, 2, . . . }, then inte- 
grating and applying an orthogonality relation to simplify the result. We note that T 
is an orthonormal set of eigenfunctions for the diffusion problem with mixed (Dirich- 
let/Neumann) and homogeneous boundary conditions, where the eigenvalues are fi n = 
(2n — 1)5. In contrast, the eigenfunctions in the problem at hand are not orthogonal 
because of the integral boundary condition at y = that leads to the more complicated 
eigenvalue equation (|5 . 19|) for which the /j, n only approach (2n — 1)5 as n — > 00. As a 
result, the eigenfunctions satisfy 



cos(^„(y- 1)) cos(pi(y - l))dy = 



I + § cos 2 (/i n ), \in = l, 
C cos(/x„) cos(^), if n ^ I. 
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If the eigenfunctions were orthogonal, then the integrals for these two cases would instead 
evaluate to \ and respectively. For the specific case with n = £ = 1, we find that 

f cos 2 ( M i(y-l))£fy« 0.4994, 
Jo 

while for n = 1 and I = 2 

/ cos(/x 1 (t/-l))cos(/x 2 (j/-l))dy« 3.701 x 10~ 4 . 
Jo 

For larger values of n and £, these integrals are even closer to the ideal values of \ 
and and therefore the eigenfunctions are very nearly orthogonal. As a result, we are 
able in practice to evaluate the series coefficients numerically by assuming that they are 
orthogonal and taking the inner solution to be 

H( + H J 1 C (y) dy ~ 2 
7(2/, t) = — + 2_^a n cos(n n (y - l))e 71 , (5.20) 

where 

a„ « 2 / (C (y) ~ 7oo)cos (^„(y - 1)) dy, (5.21) 
Jo 

and /x n are the roots of (|5.19|) . 

We remark here that a similar problem with an integral boundary condition has been 
studied by Bcilin [I], who also looked for a series solution and obtained eigenfunctions 
that are not orthogonal. However, he carried the analytical solution further by deriving 
a second set of dual eigenfunctions for an associated adjoint problem that are orthogonal 
to the original eigenfunctions. He then used both sets of eigenfunctions to calculate 
the series coefficients analytically. We have not applied Beilin's approach here because 
our problem has a more complicated integral boundary condition that leads to a time- 
dependent boundary condition in the adjoint problem for which we cannot obtain the 
eigenfunctions in the same way. 

Finally, we note that contrary to the usual approach for developing matched asymp- 
totics expansions, the inner and outer solutions in our situation involve no unspecified 
constant (s) that require matching. In particular, the inner solution for the gas concentra- 
tion tends to the constant function as t — > 00. This is also the steady state solution 
of the diffusion equation in the domain < y < 1 which coincides with the zeroth order 
term in the outer expansion as t — > 0. 

5.4 Comparison with numerical simulations 

The asymptotic solution developed in the preceding sections is now calculated using 
the same parameter values that were used in the full numerical simulations shown in 
Figures [JHU and the corresponding results are reported in Figures [7H2] respectively. In 
all cases, the inner series solution from (|5.20|) was truncated at 10 terms, while the outer 
solution is depicted for both the one- and two-term series approximations. 

Focusing first on the base case results in Figure [7t>, for very short times the inner 
concentration solution is indistinguishable to the naked eye from the computed results. 
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(a) Temperature 
(crosses - analytical solution). 



(b) Gas concentration (short time) 
(crosses - analytical solution). 
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(c) Water-ice interface. 



(d) Gas concentration at x = L. 
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Figure 7. Comparison of analytical and numerical solutions with boundary temperatures 
Z\ = T c + 0.005, T 2 = T C - 0.005, f 2 = -1. 



(a) Water-ice interface. 



(b) Gas concentration at x = L. 





Figure 8. Comparison of analytical and numerical solutions with boundary temperatures 
Zi = T c + 0.005, T 2 = T C - 0.02, f 2 = -4. 



Over longer times, the temperature and two-term series expansions for both interfacial 
position and concentration (in Figures Ek, c and d respectively) also sit directly on top of 
the computed results. The leading order concentration solution begins to deviate from the 
computed results when the boundary temperature difference is increased to T 2 — T% = — 2 
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(a) Water-ice interface. (b) Gas concentration at x — L. 





in Figure |9Jd; this reduction in accuracy derives from the fact that the Go approximation 
is constant in space, whereas the actual concentration becomes more concave in y as 
T2 — T\ increases. There is a more noticeable error in the leading order term for the 
interface position, which most evident in Figure [5£i. 

It is interesting to investigate the limitations of our asymptotic solution for more ex- 
treme values of the parameters and thereby determine under what circumstances the 
series begins to break down. For example, if the domain length is increased by several 
orders of magnitude to L = 1 cm, then there is finally a noticeable error in the two-term 
solution for concentration as shown in Figure fTOb : furthermore, the two-term asymptotic 
solution fails to adequately capture the interface position. Because it is only for such ex- 
treme values of parameters that the series approximation breaks down, we conclude that 
our approximate solutions remain accurate for the range of parameters corresponding to 
the melting of frozen sap in maple xylem cells. 
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6 Conclusions 

In this paper, we have developed a mathematical model for a three-phase free boundary 
problem that is motivated by the study of melting of frozen sap within maple trees. The 
model incorporates both melting of ice and dissolution of gas within the meltwater. We 
derive an approximate solution that captures the dynamics of the ice- water interface as a 
series expansion in the Biot number. The dissolved gas concentration exhibits variations 
over two widely disparate time scales, leading to a two-scale asymptotic solution. Com- 
parisons with numerical simulations show that the approximate solutions are accurate 
for the range of parameter values of interest in maple trees. 

There are several possible avenues for future work. First, the gas- water interfaces within 
actual xylem cells experience a large curvature, so that the interfacial surface tension will 
have a significant effect on pressure differences. We would like to include this effect, as well 
as the Gibbs-Thompson phenomenon for the variation of melting temperature across a 
curved interface which has been well-studied in the mathematical literature [16j . Secondly, 
maple trees undergo repeated daily cycles of freezing and thawing, and so the freezing 
mechanism also needs to be analysed with a daily periodic variation in the temperature. 
Finally, we would like to study further some of the technical issues surrounding the 
extension of Beilin's approach [b to the more complicated adjoint problem that derives 
from our integral boundary condition. 



Appendix A Derivation of the gas-water interface equation (|2.11[) 

Here we apply a conservation of mass argument to derive the equation (|2. 1 1|) relating 
s g w and s W i, assuming that the domain is a cylinder with constant radius r. At any time 
t, the total mass of gas is given by the integral 

/ /•«„«,(*) \ 
M g (t)=AlJ p g (s,t)ds + m w {t)\ = As gw (0)p g (0), (Al) 

while that for water is 

remit) 

M w (t) = A p w ds = 2A(s wi (t)-s gw (t))p w (A2) 

Js gw (t) 

and for ice is 

Mi(t)=A[ p i ds = 2A(L-s wi (t))p i . (A3) 

As mentioned earlier, the water and ice densities are taken to be constant. 

Since we assume that the system is closed, the sum of the three masses must be some 
constant, say .Mo, and so 

Mo =M g {t)+M w {t)+M z (t). (A 4) 

Differentiating this expression with respect to time yields 

= + A(s wi (t) - s gw (t))p w - As wl {t)pi, (A 5) 
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or simply 



,(t) 



Pi 

Pu 



i(t). 



(A 6) 



Appendix B Approximation of the eigenvalues /i J( 



Here we approximate the coefficients fi n from equation (|5 . 19[) for small values of H . Then 
the integral boundary condition reduces to a pure Dirichlet condition and we then expect 
that the eigenvalues and eigenfunctions will reduce to those of the standard separation 
of variables solution. Indeed, if = then equation (|5.19[) reduces to fi n cos(/i„) = 0, 
whose solutions are {fj,^ — (2n — l)^, n = 1,2,...}. Because we are interested in the 
case when is very small, we can make the ansatz /i„ = /i° + e„ with e„ — > 0, and 
assume further that | sin(/i„)| w 1. As a result, equation (|5.19p becomes 

H 

J' 

Using the Taylor series expansion of the cosine function centered at , (|B IP reduces to 



cos(/i„)| 



(Bl) 



2 , 



resulting in 



2 



^2 



(B2) 



(B3) 



Finally, we employ the approximation 



VT+z = 1 + - + o(z), 



to obtain 



so that 



u° 

l-'n 



IE 

1 + — 



C(m° 



\2 



H 



, 7T 2H 

fi n = (2n - lj- + 



2 C(2n-l)7r' 



(B4) 



(B5) 



Acknowledgement s 

This work was supported by a Discovery Grant from the Natural Sciences and Engineer- 
ing Research Council of Canada and a Research Grant from the North American Maple 
Syrup Council. MC was funded partially by a Fellowship from the Mprime Network of 
Centres of Excellence. 



References 



[1] S. A. Beilin. Existence of solutions for one-dimensional wave equations with nonlocal 
conditions. Electron. J. Diff. Equat., 2001(76):l-8, 2001. 



Three-phase free boundary problem with melting and dissolution 



29 



[2] H. S. Carslaw and J. C. Jaeger. Conduction of Heat in Solids. Oxford Science Publications. 

Clarendon Press, New York, second edition, 1988. 
[3] M. Ceseri and J. M. Stockie. A mathematical model for sap exudation in maple trees 

governed by ice melting, gas dissolution and osmosis. To appear in SIAM J. Appl. 

Math., 2012. 

[4] J. Crank. The Mathematics of Diffusion. Clarendon Press, 1956. 
[5] J. Crank. Free and Moving Boundary Problems. Clarendon Press, New York, 1984. 
[6] A. Friedman. Free boundary problems for parabolic equations. I. Melting of solids. J. Math. 
Mech., 8:499-517, 1959. 

[7] A. Friedman. Free boundary problems for parabolic equations. III. Dissolution of a gas 

bubble in liquid. J. Math. Mech., 9:327-345, 1960. 
[8] A. Friedman. Variational Principles and Free-Boundary Problems. John Wiley & Sons, 

New York, 1982. 

[9] A. Friedman. Free boundary problems in science and technology. AMS Notices, 47(8) :854- 
861, 2000. 

[10] R. M. Furzeland. A comparative study of numerical methods for moving boundary prob- 
lems. J. Inst. Math. Appl, 26(4) :41 1-429, 1980. 

[11] S. C. Gupta. The Classican Stefan Problem: Basic Concepts, Modelling and Analysis, 
volume 45 of North-Holland Series in Applied Mathematics and Mechanics. Elsevier, 
Amsterdam, 2003. 

[12] W. Huang and R. D. Russell. Adaptive Moving Mesh Methods, volume 174 of Applied 

Mathematical Sciences. Springer, New York, 2011. 
[13] P. Huyakorn, S. Panday, and Y. Wu. A three-dimensional multiphase flow model for assesing 

NAPL contamination in porous and fractured media, 1. Formulation. J. Contam. Hydrol., 

16(2): 109-130, 1994. 

[14] J. B. Keller. Growth and decay of gas bubbles in liquids. In Proceedings of the Symposium 

on Cavitation in Real Liquids (General Motors Research Laboratory, Warren, MI), pages 

19-29. Elsevier, New York, 1964. 
[15] W. Konrad and A. Roth-Nebelsick. The dynamics of gas bubbles in conduits of vascular 

plants and implications for embolism repair. J. Theor. Bio., 224:43-61, 2003. 
[16] S. Luckhaus. Solutions for the two-phase Stefan problem with the Gibbs-Thomson Law for 

the melting temperature. Euro. J. Appl. Math., 1:101-111, 1990. 
[17] J. Milburn and P. O'Malley. Freeze-induced fluctuations in xylem sap pressure in Acer 

pseudoplatanus: A possible mechanism. Can. J. Bot., 62:2100-2106, 1984. 
[18] M. S. Plesset and A. Prosperetti. Bubble dynamics and cavitation. Annu. Rev. Fluid Mech., 

9:145-185, 1977. 

[19] K. Singh and R. K. Niven. Non-aqueous phase liquid spills in freezing and thawing soils: 

Critical analysis of pore-scale processes. Crit. Rev. Environ. Sci. Technol, 2012. In press, 

DOI: 10.1080/10643389.2011.604264. 
[20] L. N. Tao. On solidification problems including the density jump at the moving boundary. 

Quart. J. Mech. Appl. Math., 32(2): 175-185, 1979. 
[21] G. G. Tsypkin. Mathematical models of gas hydrates dissociation in porous media. Ann. 

New York Acad. Set., 912(l):428-436, 2000. 
[22] M. T. Tyree and J. S. Sperry. Vulnerability of xylem to cavitation and embolism. Annu. 

Rev. Plant Physiol. Plant Mol. Biol, 40:19-38, 1989. 
[23] D. G. Wilson. One dimensional multi-phase moving boundary problems with phases of 

different densities. Technical Report CSD-93, Oak Ridge National Laboratory, January 

1982. 

[24] W. Xu. Modeling dynamic marine gas hydrate systems. Amer. Mineral, 89(8-9) :1271-1279, 
2004. 



